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Abstract 

Spatial adaption procedures for the accurate and effi- 
cient solution of steady and unsteady inviscid flow prob- 
lems are described. The adaption procedures were developed 
and implemented within a two-dimensional unstructured- 
grid upwind-type Euler code. These procedures involve 
mesh enrichment and mesh coarsening to cither add points 
in high gradient regions of the flow or remove points where 
they are not needed, respectively, to produce solutions of 
high spatial accuracy at minimal computational cost The 
paper gives a detailed description of the enrichment and 
coarsening procedures and presents comparisons with alter- 
native results and experimental data to provide an assessment 
of the accuracy and efficiency of the capability. Steady and 
unsteady transonic results, obtained using spatial adaption 
for the NACA 0012 airfoil, are shown to be of high spatial 
accuracy, primarily in that the shock waves are very sharply 
captured. The results were obtained with a computational 
savings of a factor of approximately fifty-three for a steady 
case and as much as twenty-five for the unsteady cases. 

Introduction 

Considerable progress has been made over the past 
decade in developing methods of dynamically-adapting com- 
putational meshes based on the numerical solution of partial 
differential equations. 1 These methods are being developed 
to produce higher spatial accuracy in such solutions more ef- 
ficiently. Spatial accuracy is obviously important when mod- 
elling continuous equations with a discrete set of points. It 
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is generally understood, of course, that accuracy is improved 
when the number of mesh points in a fixed computational 
domain is increased. Associated with an increase in the 
number of mesh points, however, are increased computer 
run times and memory costs, which are proportional to the 
number of mesh points. Hence, for efficiency, it is important 
to enrich meshes locally based on the numerical solution, in 
contrast to using globally fine meshes, to minimize the total 
number of mesh points and hence minimize the cost for a 
given spatial accuracy. The methods of mesh adaption can 
be separated into three general categories: (1) mesh regen- 
eration, (2) mesh movement, and (3) mesh enrichment. The 
first method, mesh regeneration, places the work of adapt- 
ing the mesh on the mesh generation program rather than 
on the actual numerical solution procedure of the governing 
equations. In this method, a solution is first obtained, and 
regions of relatively large discretization errors are detected. 
A new mesh is then generated to concentrate points in re- 
gions where the large discretization errors occur. This new 
mesh may contain more or fewer points than the original 
mesh. For the second method, mesh movement, the number 
of points in the computational domain remains fixed. To re- 
solve more accurately the solution spatially, these points are 
moved into regions where solution gradients are relatively 
large. In general, this can be accomplished by two different 
ways. The first way models the mesh as a spring network, 
where points are joined by linear springs with spring stiff- 
nesses proportional to solution gradients. The mesh is then 
allowed to move into the relatively high gradient regions to 
produce effectively a locally finer mesh. The second way 
uses forcing functions in aFoisson-equation grid generator to 
redistribute points. Either method of mesh movement is eas- 
ily implemented within existing solution algorithms because 
only the locations of the existing mesh points are changed. 
The final method of spatial adaption is mesh enrichment. In 
this method points arc added to regions of relatively large so- 
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lution error by dividing locally the cells which make up the 
mesh or by embedding finer meshes in these regions. This 
method differs from mesh regeneration and movement in that 
the mesh is made finer in local regions while the global mesh 
topology remains the same. The method of mesh enrichment 
is also generally regarded as having advantages over regen- 
eration and movement, especially for transient problems, be- 
cause of the higher degree of flexibility in being able to add 
points where they are needed and to remove points where 
they are not needed. The disadvantage, however, is that 
the implementation of mesh enrichment involves significant 
modification to existing solution algorithms. 

For the Euler and Navier-Stokes equations, computa- 
tional fluid dynamics algorithms arc being developed based 
on spatial adaption methods. With these equations, rel- 
atively large spatial discretization errors may be encoun- 
tered with flow features such as shock waves, shear layers, 
boundary layers, and expansion fans. These flow features 
can be resolved more accurately using the adaption meth- 
ods mentioned above, Nakahashi el al, 2 for example, has 
used tension and torsional springs to move the mesh into 
regions where relatively large spatial discretization errors 
occur. This mesh movement approach showed consider- 
able versatility for the problems treated. However, various 
constants were needed to control orthogonality and smooth- 
ness, and direct control of an optimal mesh adaption proce- 
dure was generally not possible. Further examples of spatial 
adaption methods include the work of Usab and Murman. 3 In 
Ref. 3, embedded meshes of quadrilateral cells and nodes 
were used in regions of the mesh where shock waves oc- 
curred, This approach improved the spatial accuracy of the 
numerical method which resulted in highly accurate solu- 
tions for steady flow problems. Dannenhoffer and Baron 4,5 
extended the work in this area using irregularly shaped em- 
bedded regions, which were coupled to the base mesh by 
a multiple-grid solution algorithm. Several other examples 
of spatial adaption include methods which use flow solvers 
based on unstructured triangular and tetrahedral meshes in 
two and three dimensions, respectively. Pcraire et al. 6,7 used 
mesh regeneration coupled with a finite element solution al- 
gorithm to sharply capture shock waves and complex shock 
structures. Lohner et al. 8 ’ 13 have developed a procedure 
to locally enrich the mesh for transient flow problems by 
dividing elements which make up a base mesh to capture 
shock waves. Further, in this procedure, elements may be 
removed (coarsened) from the mesh if they are not necessary 
to produce a given level of spatial accuracy. 

With respect to solution algorithms based on unstruc- 
tured meshes, the results published by the authors 14 demon- 
strated that these algorithms produce steady and unsteady 
solutions of comparable accuracy to results obtained using 
structured-grid solution algorithms. The two sets of results 
presented in Ref. 14 were obtained using meshes of compa- 
rable density, and mesh adaption was not used. Solutions of 
higher spatial accuracy are indeed possible through the use 


of mesh adaption. Therefore, the purpose of this paper is to 
report on modifications to the two-dimensional unstructured- 
grid upwind-type Euler code of Batina 15 to include mesh en- 
richment and coarsening procedures. The objectives of the 
research arc as follows: (1) to develop time -accurate enrich- 
ment and coarsening procedures for spatial adaption, (2) to 
test the procedures by performing steady and unsteady calcu- 
lations for a variety of cases, (3) to determine the accuracy of 
the spatially adapted solutions by making comparisons with 
published solutions produced by alternative methods and ex- 
isting experimental data and (4) to assess the efficiency of 
the spatially adapted solutions by making comparisons of re- 
quired computer resources. The eventual goal is to develop 
a highly accurate and efficient solution algorithm for the Eu- 
ler and Navier-Stokes equations for aeroelastic analysis of 
complex aircraft configurations. The paper gives a detailed 
description of the mesh enrichment and coarsening proce- 
dures and gives a brief description of the solution algorithm 
of Ref. 15 for completeness. Steady and unsteady transonic 
results are presented for the NACA 0012 airfoil to demon- 
strate applications of the spatial adaption procedures. The 
unsteady flow results for the NACA 0012 airfoil were ob- 
tained for the airfoil pitching harmonically about the quartei 
chord. 


Upwind-Type Euler Solution Algorithm 

The unsteady Euler equations arc solved using flic two- 
dimensional upwind-type solution algorithm developed by 
Batina. 15 The solution algorithm, which is a cell-centered 
finite-volume scheme, 1517 uses upwind differencing based 
on flux-vector splitting similar to upwind schemes developed 
for use on structured meshes. The present unstructured grid 
algorithm is thus referred to as an upwind-type solution al- 
gorithm. The spatial discretization of this algorithm involves 
a so-called flux-split approach based on the flux-vector split- 
ting of van Leer. 18 The flux-split discretization accounts for 
the local wave-propagation charactcristics of the flow and 
captures shock waves sharply with at most one grid point 
within the shock structure. A further advantage is that the 
discretization is naturally dissipative and consequently does 
not require additional artificial dissipation terms or the ad- 
justment of free parameters to control the dissipation. How- 
ever, in calculations involving higher-order upwind schemes, 
oscillations in the solution near shock waves are expected to 
occur. To eliminate these oscillations, flux limiting is usually 
required. In the present study, a continuously differe ntia ble 
flux limiter was employed. 15-17 

The Euler equations are integrated in time using an 
implicit time- integration scheme involving a Gauss-Seidel 
relaxation procedure. 15 The relaxation procedure is imple- 
mented by reordering the elements that make up the unstruc- 
tured mesh from upstream to downstream. The solution is 
obtained by sweeping two times through the mesh as dictated 
by stability considerations. The first sweep is performed in 
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the direction from upstream to downstream and the second 
sweep is from downstream to upstream. For purely super- 
sonic flows the second sweep is unnecessary. This relaxation 
scheme is stable for large time steps and allows the selection 
of the step size based on the temporal accuracy required for 
the problem being considered, rather than on the numerical 
stability of the algorithm. Consequently, very large time 
steps may be used for rapid convergence to steady state and 
an appropriate step size may be selected based on temporal 
convergence for unsteady cases, independent of numerical 
stability issues. 

Weighted Averaging 


using the surrounding nodal values fa and the weights w,. 
The weights were derived by first noting that the Laplacian 
of a linear function is exactly zero. Since this is a desirable 
property to require of the computed Laplacian, the weights 
w, in Eq. (2) were derived with this objective. Therefore, 
the spatial locations of <j>, were substituted into Eq. 

(2) resulting in 

n 

L(*o) = “'•(*• - *o) = o (3) 

1 = 1 

n 

L(y o) = w '( y ' - M>) = o (4) 

1=1 


A weighted averaging procedure is used in the two- 
dimensional upwind-type Euler solution algorithm to inter- 
polate the cell-centered values of the primitive variables to 
the nodes. As illustrated in Fig. 1, the weighted average of 
the cell-centered values q % surrounding a node is given by, 

n 

E »"< 0» 

m (1) 

£ m 

»=i 

where the nodal value q 0 is computed using the cell weights 
Wi. Several different weights have been used in this proce- 



Fig. 1 Difference star for computing the weighted average 
q 0 of surrounding cell-centered values q x ... n , 

dure including the areas of surrounding cells and the recip- 
rocal distance from the node to the locations of cell centers. 
Because of the disparity in cell sizes that occur when using 
spatial adaption, the most accurate weighted averaging is 
desired. Therefore, a different approach for computing the 
weights has been developed based on the work of Holmes 
and Connell. 19 In Ref. 19, artificial dissipation for a node- 
based Navier-Stokes solution algorithm was computed using 
a Laplacian operator. The Laplacian L(<f> 0 ) at a node was 
computed by 


The weights were then determined by defining 


rr 1 + (5) 

where 

n 

<" = £( A ^) 2 ( 6 ) 

1=1 

is a cost function. The cost function C was minimized to 
keep the weights close to unity by solving an optimization 
problem given the constraints of Eqs. (3) and (4). The opti- 
mization problem was solved using the method of Lagrange 
multipliers where A u>, was given by 

Au>i = A X (xi - x 0 ) + A y(y { - y 0 ) (7) 


The solution yielded Lagrange multipliers, defined by 


where 
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Note that the weights were computed solely from the dis- 
(2) tancc between node locations and no special treatment was 
required near mesh boundaries. 
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Modifications were made to this procedure for comput- 
ing the Laplacian in order to develop a new procedure for 
computing the weights in Eq. (1). The first modification 
was to express the Laplacian at the nodes using surrounding 
cell -centered values, as illustrated in Fig. 1, rather than the 
surrounding nodal values. Now the Laplacian at the node 
can be expressed by 

n 

L (<lo) = “><(?< - ?o) (15) 

i — 1 

where q 0 represents the unknown value at the node and q x 
are the known values at the surrounding cell-centers. The 
weights Wi in Eq. (15) are now computed using the node 
location (x 0x yo) and the surrounding cell-center locations 
(**!»)• Recalling that the computed Laplacian for a linear 
function is zero, then for a linear variation of q 

n 

w <(n -?o) = o (16) 

» = 1 

which can be rewritten in the form given by Eq. (1). 
Therefore, using the weights from the Laplacian procedure, 
Eq. (1) is an exact inteipolation of the nodal value q 0 for a 
linear variation of cell-centered values q t . For a nonlinear 
variation of cell-centered values q t , the weighted averaging 
procedure using the weights from the Laplacian procedure 
is second order accurate in space. Although the procedure 
is computationally more expensive than area or reciprocal 
distance weighted averaging, the increased cost is relatively 
small compared to the cost of solving the Euler equations. 

Spatial Adaption Procedures 

In this section, the spatial adaption procedures are de- 
scribed in detail. These descriptions include explanations of 
the flow feature detection and the procedures used to enrich 
and coarsen the mesh. 

Flow Feature Detection 

The first step of the spatial adaption procedure is the 
detection of regions of relatively large discretization error 
so that the computational mesh can be locally enriched to 
improve the spatial accuracy or locally coarsened to reduce 
the computational costs. These regions generally occur near 
flow features such as shock waves, stagnation points, slip 
lines, and expansion fans for the Euler equations, where the 
dominant flow feature for the cases considered in this study 
is a shock wave. 

There are a number of flow parameters that can be 
used for enrichment indicators based on the detection of 
shock waves. Parameters such as density, pressure, or total 
velocity are useful since these quantities are discontinuous 
through shocks. For example, first or second, divided or 


undivided differences in one of these parameters, similar 
to the work by Dannenhoffer and Baron, 4 can be used to 
delect shock waves. A specific example of an enrichment 
indicator is the magnitude of the density gradient | Vp| which 
is often used to detect shock waves for steady flow problems. 
However, for unsteady flows, a measure of the temporal as 
well as the spa tial variat ion of the solution is needed. For 
this reason, the absolute value of the substantial derivative 
of density |^| was used as an enrichment indicator, since 

it is the sum of the local rate of change as well as the 
convective rate of change u Vp of density. The substantial 
derivative of density can also be written as 

Dp dp _ 

S = 5 4V "-' vi < 17) 

where the first two terms on the right-hand-side of Eq. (17) 
are simply the continuity equation. The third term can 
be viewed as a simplified continuity equation that assumes 
steady, incompressible, linear flow. Therefore, the differ- 
ence between these two continuity equations results in a 
measure of the 1) unsteadiness, 2) compressibility, and 3) 
nonlinearities of the flow. Since the first two terms on the 
right-hand-side of Eq. (17) are zero from the continuity 
equation, the substantial derivative of density can be rewrit- 
ten as 

-£ = -pV u (18) 

For unsteady flows, the substantial derivative works well for 
detecting developing shock waves, whereas the more com- 
monly used enrichment indicators based on the instantaneous 
solution such as first or second differences in density tend 
to miss the initial shock wave formation. This is especially 
true for cases where the shock waves periodically appear and 
disappear in time. Other indicators, based on the substantial 
derivative of x-momentum pu, y-momentum pv, and energy 
e have also been investigated. These substantial derivatives 
are written as 

Dpu _ j dpu d(pu 2 + p) puv \ 

Dt ~ \ ~dt~ + di + ~dy ' 

Dpv _ f dpv dpuv d(pv 2 4- p) 

Dt | dt + dx + dy 

De _ fde. ^ d(e + p)u d(e + p)v ) 

Dt \ dt dx dy / ( 21 ) 

— {V ( pit) + eV «} 

where once again each derivative results in the difference be- 
tween a complete conservation equation and an approximate 
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equation. The first expression in brackets on the right-hand- 
side of Eqs. (19)-(21) is zero for the Euler x-momentum, 
y-momentum, and eneigy equations, respectively. The last 
expression in brackets is an approximate equation of the re- 
spective conservation laws. For example, the second equa- 
tion in brackets in Eq. (19) can be derived from the conser- 
vation of x-momentum equation if the equation is simplified 
by assuming the local rate of change and the convective 
rate of change u V of x-momentum is zero. Therefore, 
as one might expect, if the absolute value of Eq. (19) is 
used as an enrichment indicator it would delect regions of 
the flow where the x-momentum is most rapidly changing. 
Similarly, Eqs. (20) and (21) would detect the rale of change 
of y-momentum and energy, respectively. To evaluate each 
expression as an enrichment indicator, the substantial deriva- 
tives of the conserved variables were computed for a tran- 
sient shock problem. In this study, the performance of each 
of the substantial derivatives in detecting shock waves was 
found to be similar. It is not surprising, however, since all 
of the indicators involve the same term V u. 

To detect viscous flow features, the idea of using 
the difference between complete equations and approxi- 
mate equations can be expanded upon by considering the 
differences between the momentum or energy equation of 
the Navier-Stokes equations and an approximate inviscid 
equation. For example, using the momentum equation of 
the Navier-Stokes equations, the substantial derivative of x- 
momentum can be written as 


Dpu _ f dpu d (pu 2 + P - T rx ) d(puv - r xy ) 
Dt ~ | dt + dx + ~dy 



dp _ 

+ puV ■ u 
ox 


OTxi 

dx 


d±ry \ 

dy J 


( 22 ) 

This, of course, is similar to Eq. (19) except that the 
viscous fluxes ( ^ ) have been added and subtracted 
on the right-hand-side. Therefore, using the ideas above, 
an enrichment indicator E for detecting inviscid as well 
as viscous flow features can be constructed by eliminating 
the viscous fiuxcs from the second expression enclosed in 
brackets of Eq. (22) and is written as 


E = 


Op J-7 

-f puV • u 
Ox 


(23) 


Although the enrichment indicator in Eq. (23) is not used 
in this study, it is mentioned to demonstrate a procedure 
for constructing an indicator based on the governing viscous 
flow equations. The enrichment indicator used in this study, 
however, was the absolute value of the substantial derivative 
of density (Eq. (18)) since the Euler equations were used. 


Mesh Enrichment 

Generally, mesh enrichment is performed by starting 
with a relatively coarse mesh of cells and then subdividing 


these cells until a given level of spatial accuracy has been 
obtained. To prevent cells from being enriched too many 
times near flow discontinuities such as shock waves, an up- 
per bound is placed on the number of times a cell can be 
divided. For transient problems the mesh enrichment proce- 
dure may be performed at each time step of the integration 
of the governing flow equations or it may be performed once 
every set number of time steps. 



(b) lypc-2 clement enrichment. 



(c) restriction on further enrichment of type-2 elements. 



(d) details on further enrichment of type-2 elements. 


Fig. 2 Diagrams illustrating mesh enrichment possibilities. 


Mesh enrichment is performed by using the enrichment 
indicator to determine if a cell is to be subdivided into 
smaller cells. To accomplish this, the enrichment indicator is 
computed for each cell and compared with a preset threshold 
value to determine whether a cell should be subdivided. If 
the threshold is exceeded, a new node is created at the 
midpoint of each edge of the triangular cell, and the cell 
is subdivided into four smaller cells. Special care must be 
taken, however, when an edge that is to be bisected lies 
on a boundary of the mesh, since the midpoint of the edge 
docs not generally lie on the boundary. In this case, the 
location of the new node is determined by using a parametric 
spline of the boundary coordinates. Further, the values of 
the flow variables for the new cells are determined from 
a linear interpolation of conserved variables located at the 
nodes and the cell centers of the original cells. 

For a given cell to be enriched, either one edge or all 
three edges arc bisected. In the event that only two edges 
arc marked to be bisected, the third edge is automatically 
bisected to prevent the creation of highly skewed or stretched 
cells. Each time the mesh is enriched, a cell may be divided 
in one of two ways, as shown in Figs. 2(a) and 2(b). The 
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first way, shown in Fig. 2(a), results when all three edges 
of a cell have been marked for enrichment. In this situation 
the cell is divided into four new cells where the vertices of 
the inner cell are in general midpoints of edges that make 
up the original cell. The original cell is thus referred to as a 
type-4 element since after enrichment it becomes four new 
triangular cells. The second way, shown in Fig. 2(b), occurs 
if only one edge of the base cell is marked for enrichment. 
In this situation, the marked edge is bisected, and two new 
cells are formed. The original cell is thus referred to as a 
type-2 element since after enrichment it becomes two new 
cells. New cells formed from a type-4 element may be 
enriched further. However, to prevent highly stretched cells, 
cells formed from a type-2 element are restricted from being 
divided further as indicated in Figs. 2(c) and 2(d). For cells 
from a type-2 element, if any of the five edges that make up 
the two new cells are marked for enrichment, the original 
cell is made into a type-4 element as shown in Fig. 2(c). If 
in addition to this, either or both of the bottom two edges 
are marked (the lower left, right, or both), cells of the type-4 
element are further enriched accordingly, as shown in Fig. 
2(d). 

Mesh Coarsening 

Generally, mesh coarsening removes added nodes and 
cells from previously enriched meshes to delete them from 
local regions of the mesh where certain flow features are 
no longer present. This procedure is necessary to efficiently 
adapt meshes to the numerical solution of the governing 
flow equations in order to minimize computational cost. 
The mesh coarsening procedure starts by marking all of 
the cells that do not have edges that are marked to be 
bisected from the mesh enrichment procedure. Figure 3 
shows an example of the coarsening procedure where the 
dashed lines represent cells formed by a previous enrichment 
of the mesh. In Fig. 3(a), cells which are candidates for 
removal are denoted by a 'T\ The marked cells are then 
used to determine nodes that are candidates for removal. 
As shown in Fig. 3(b), the candidate nodes that may be 
removed are nodes that are surrounded by cells that arc 
candidates for removal (identified by the ’’ones” in Fig. 3(a)). 
Cells that are to be removed from the list of candidates are 
determined by searching the tree data structure that stores 
the mesh enrichment history. Cells for removal are cells that 
came from type-2 or type-4 elements and are subsequently 
marked for removal. Once the nodes and elements have been 
selected they are subjected to a final simultaneous evaluation 
for actual removal from the mesh. Each time the mesh 
is coarsened, cells and nodes may be removed in one of 
several ways as shown in Fig. 4. (In the figure, nodes to be 
removed are indicated by the open circles and fixed nodes 
are indicated by the closed circles.) For a type-4 element, 
if all three of the nodes that form the inner triangle arc 
candidate nodes, then the nodes arc removed eliminating the 
inner triangle as shown in Fig. 4(a). This leaves only the one 


original cell that was previously divided into four. Similarly, 
if two of the three nodes that form the inner triangle are 
candidate nodes, the two nodes are removed and a type-2 
element is formed as shown in Fig. 4(b). If only one of the 
candidate nodes is marked then nothing is done. For a type- 
2 element there is only one node that may be removed which 
is the midpoint of a previously bisected edge. Removal of 
this node leaves only the cell that was originally divided 
into two as shown in Fig. 4(c). 


0 fixed cell • fixed node 

1 candidate cell to be O candidate node to be 

removed removed 




(a) identify cells to (b)identify nodes to (c) final mesh, 
be removed. be removed . 

Fig. 3 Diagrams illustrating details of the mesh 
coarsening procedure. 




(b) two nodes removed from type-4 element 



Fig. 4 Diagrams illustrating mesh coarsening possibilities. 
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Temporal Convergence Considerations 

Careful consideration must be given to the time-step 
chosen to ensure temporal convergence of unsteady prob- 
lems, especially when spatial adaption procedures are used. 
In general the temporal convergence of the numerical solu- 
tion of a time accurate problem is difficult to quantify and 
is related to the nondimensional time-step A t 9 a characteris- 
tic length (A/s), and the local maximum characteristic wave 
speed (|ti| -f a) of each cell in the computational mesh. This 
relation is given by the Courant-Friedrichs-Lewy ( CFL ) 
number for each cell, 


from Eq. (24), the lime-step is proportional to the charac- 
teristic length (A/s), the allowable time-step is halved when 
the cell is enriched one level. 


At <x 


A/4 

s/2 


1 A 

2 8 


(25) 


For unsteady problems involving a pitching airfoil in forced 
harmonic motion, the effect of enrichment on the number 
of time-steps per cycle of motion can also be demonstrated. 
First the number of time-steps per cycle of motion N is 
related to the lime-step by 


CFL sa A<-^-(|«| + a) (24) 

where s is the length of a cell edge, A is the area of the 
cell, u is a local flow speed normal to a cell edge, and a 
is the local speed of sound. For unsteady problems, expe- 
rience has shown that the local CFL number must remain 
below approximately 20 for solutions to be converged in 
time. However, the local allowable CFL number is highly 
dependent on the time -marching procedure used. For exam- 
ple, either explicit or implicit timc-marching procedures may 
be used to integrate the governing flow equations in time for 
unsteady problems. Since the allowable CFL number based 
on a stability analysis of an explicit timc-marching procedure 
is usually small compared to 20, the time-step is restricted to 
be relatively small. Therefore, the numerical stability of an 
explicit time-marching procedure dictates the time-step that 
can be used. For an implicit time-marching procedure, how- 
ever, numerical stability of the procedure is, in general, not 
an issue. Hence, the time-step should be based on minimiz- 
ing the effects of linearization errors, minimizing the effects 
of relaxation errors (or factorization errors in the case of 
a factored scheme), and obtaining a temporally converged 
solution. Irrespective of the time-marching procedure used, 
the spatial adaption procedures have a significant effect on 
the allowable lime-step when the CFL number is restricted 
by either the temporal accuracy of the solution or the nu- 
merical stability of the solution algorithm. This effect can 
be demonstrated in a simple example of mesh enrichment of 
a single cell. Suppose, for example, the local CFL number 
for a cell is restricted to 20 for a temporally converged so- 
lution and the cell is enriched one level, where it is divided 
into four smaller cells each with area A/4 as shown in Fig. 
5. Likewise, each edge of the cell is bisected resulting in 
each of the new cells having edge lengths of s/2. Since, 




Fig. 5 Diagram illustrating change in characteristic length 
and area of an enriched computational cell. 


N = 


7 T 

k Moo At 


(26) 


where k is the reduced frequency and M is the freestream 
Mach number (since At is nondimensionalized by the 
freestream speed of sound). Therefore, the number of time- 
steps per cycle of motion is inversely related to the time- 
step. 


N cx 


a7 


(27) 


Since the time-step is halved as the mesh is enriched one 
level, the number of lime-steps required per cycle of motion 
doubles for a similar level of temporal convergence. Hence, 
the relationship between the number of time-steps per cycle 
of motion N required on a mesh to the number of time-steps 
per cycle of motion N n on the same mesh using n levels of 
mesh enrichment can be stated as 


N u k2”N (28) 

As an example, suppose an unsteady problem on a coarse 
mesh requires 1250 time-steps per cycle of motion for a 
temporally converged solution. The same mesh and case us- 
ing three levels of enrichment would require approximately 
10000 time-steps per cycle of motion. 


Results and Discussion 

Calculations were performed for the NACA 0012 airfoil 
to assess the accuracy and efficiency of the spatial adaption 
procedures. Steady and unsteady airfoil results are presented 
for comparison with published solutions using alternative 
methods and experimental data to determine the accuracy of 
the results. Timing comparisons are also made between re- 
sults obtained using globally (estimated) and locally adapted 
meshes to determine the computational savings obtained by 
using the spatial adaption procedures. 

Steady Results 

Steady flow results were obtained for the NACA 0012 
airfoil for comparison with results reported by Pulliam and 
Barton 20 and several other alternative methods for an in vis- 
cid case chosen by the AGARD Working Group 07, a sub- 
panel of the AGARD Fluid Dynamics Panel. The case of 
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Fig. 7 Comparison of pressure contour lines (Ap = 0.025) 
0012 airfoil at = 0.8 and a 0 = 1.25°. 

interest is AGARD01, a NACA 0012 airfoil at a freestream 
Mach number of = 0,8 and an angle of attack of a 0 
= 1.25°. For this case, the starting mesh and final spa- 
tially adapted mesh are presented along with corresponding 
pressure contours for comparison with results of Ref. 20. 
The starting coarse mesh for the calculation, shown in Fig. 
6(a), contains 1854 nodes and 3628 cells and extends 20 
chordlengths to a circular outer boundary. The solution 
was obtained by adapting the mesh to the solution every 
500 iterations for the first 2000 iterations (four levels of en- 



from solutions obtained for the NACA 


richmcnt), and then the solution was converged to machine 
zero on ihe final adapted mesh. The final mesh after four 
levels of enrichment, shown in Fig. 6(b), contains 8947 
nodes and 17653 cells. A comparison of pressure contours 
(A p - 0.025) from the solutions obtained on the original 
mesh and the adapted mesh are shown in Figs. 7(a) and 
7(b), respectively. For this case, there is a relatively strong 
shock wave on the upper surface of the airfoil near 62% 
chord and a relatively weak shock wave on the lower sur- 
face near 30% chord. Ihc comparison between the two scls 
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(a) mesh. 



(b) pressure contours. 


Fig. 8 Results from Pulliam and Barton 20 for the NACA 0012 airfoil at M IXJ — 0.8 and < *o — 1-25 . 


of contour lines reveals a considerable improvement in the 
shock resolution when spatial adaption is used. For further 
comparison, the mesh and pressure contour lines from the 
solution of Ref. 20 are shown in Fig. 8, which is believed 
to be one of the most accurate of the results generated by the 
alternative methods. Results from Ref. 20 were computed 
on a structured mesh of C-typc topology containing 560 x 
65 mesh points as shown in Fig. 8(a). This mesh indicates 
that points were clustered above and below the airfoil in the 
shock regions to produce a more accurate result. A compari- 
son of the pressure contour lines from Pulliam and Barton, 20 
shown in Fig. 8(b), with the contour lines shown in Fig. 
7(b), obtained with the present spatial adaption procedure, 
indicates excellent agreement which verifies the accuracy of 
the adaption procedures. Finally, comparisons of computed 
lift, drag, and moment coefficients from the structured mesh 
solutions reported in Ref. 20 and the present results arc 
given in Tbble 1. The present results are shown at the 
bottom of Table 1 and the coefficients listed from solutions 
5 and 5* are from Pulliam and Barton. 20 The remaining so- 
lutions are from alternative methods and are included for 
completeness. Solution 5* is an improved result in compar- 
ison with solution 5, where improved boundary conditions 
and a finer mesh were used. A comparison of coefficients 
from the present solution with those of solution 5*. indicates 
agreement to within 1% for the lift and drag coefficients and 
to within 5% for the moment coefficient. 

The computational savings for the steady flow prob- 
lem can be estimated by comparing the number of cells 
in the adapted mesh to the number of cells in a globally 
fine mesh. The total number of cells for the globally fine 
mesh is computed by multiplying the number of cells in 
the starting coarse mesh by four (since each cell is divided 


Table 1 Comparisons of computed lift, drag, and moment 
coefficients from the structured mesh solutions 
reported in Ref. 20 and the present results 
obtained using spatial adaption for the NACA 
0012 airfoil at Moo = 0.8 and « 0 = 1.25°. 


Soln# 

Type 

Mesh Size 

Cl 

Cd 

Cm 

1 

c 

189X25 

0.3642 

0.0225 

-0.0376 

2 

c 

188X24 

0.3736 

0.0244 

-0.0411 

3 

o 

158 X 23 

0.3463 

0.0223 

-0.0358 

5 

c 

249 X 67 

0.3661 

0.0229 

-0.0430 

6 

o 

192 X 39 

0.3474 

0.0221 

-0.0374 

8 

o 

128 X 28 

0.3500 

0.0221 

-0.0370 

9 

o 

320 X 64 

0.3632 

0.0230 

-0.0397 

5* 

c 

560 X 65 

0.3618 

0.0236 

-0.0411 

present 

u 

8947N, 17653C 

0.3600 

0.0238 

-0.0390 


into four smaller cells) raised to a power represented by 
the number of enrichment levels. Therefore, a globally fine 
mesh would contain 928768 cells (3628 x 4 1 ), which in 
comparison with the final adapted mesh containing 17653 
cells, results in a computational savings of approximately 
fifty-three. If timing comparisons arc made with solutions 
obtained using structured grid algorithms, the computational 
overhead of the indirect addressing of the unstructured grid 
algorithm must be accounted for. For example, the compu- 
tational work for steady-state solutions obtained using un- 
structured grid algorithms have been shown to be 2 to 5 
limes more expensive than those obtained using structured 
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a(T) = 3.49°, kt = 214° 


a(t) = 2.43°, kt - 265° 


a(t) - 2.67°, kt = 296° 


a(t) = 4.28°, kt = 346° 


Fig. 9 Instantaneous meshes produced by the spatial adaption procedure for the NACA 0012 airfoil pitching 
harmonically at = 0.599, a 0 = 4.86°, a x = 2.44°, and k = 0.0814. 



Fig. 10 Instantaneous contour lines (Ap = 0.02) from the spatially adapted solution for the NACA 0012 airfoil 
pitching harmonically at = 0.599, « 0 = 4.86°, a, = 2.44°, and k = 0.0814. 
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grid algorithms. 21 Nonetheless, when spatial adaption pro- 
cedures are used with unstructured grid algorithms, an order 
of magnitude computational savings can be achieved over 
solutions obtained using structured grid algorithms for com- 
parable spatial accuracy. 

Unsteady Results 

Unsteady results are presented for the NACA 0012 air- 
foil using spatial adaption for AGARD cases 3 and 5, pro- 
posed by the AGARD Structures and Materials Panel. 22 Both 
cases involve the NACA 0012 airfoil pitching harmonically 
about the quarter chord at a reduced frequency based on 
semichord of k = 0.0814. The calculations were performed 
for three cycles of motion to obtain a periodic solution us- 
ing 10000 time-steps per cycle starting with the same coarse 
mesh (Fig. 6(a)) that was used for the non-adapted steady 
flow result During the course of the calculation, the mesh 
was spatially adapted every 10 time-steps using three levels 
of enrichment 

AGARD Case 3 

Results were obtained for the airfoil pitching with an 
amplitude of a\ « 2.44° at = 0.599 and a 0 = 4.86° 
(referred to as AGARD case 3). Figure 9 shows the in- 
stantaneous adapted meshes and Fig. 10 shows the corre- 
sponding instantaneous density contour lines (A p = 0.02). 


The instantaneous meshes and density contour lines during 
the third cycle of motion were plotted at eight points in 
time. In each plot, the instantaneous pitch angle a (r) and 
the instantaneous angular position in the cycle kr are noted. 
The instantaneous meshes (Fig. 9) clearly indicate the en- 
richment in regions near the shock wave on the upper sur- 
face of the airfoil and near the stagnation points. Figure 9 
also shows that regions of the mesh near the moving shock 
and in the wake of the airfoil have been enriched and then 
coarsened throughout the cycle of motion. Specifically, the 
instantaneous mesh in Fig. 9 at kr = 26° shows that the 
enrichment indicator works well for detecting the develop- 
ing shock wave on the upper surface as is further shown 
in the density contours (Fig. 10). These density contours 
during the cycle demonstrate the ability of the spatial adap- 
tion procedures to produce sharp transient shock waves. The 
corresponding surface pressure distributions during the third 
cycle of motion arc shown in Fig. 1 1 for comparison with 
the experimental data of Ref. 23. In each pressure plot 
the instantaneous pitch angle a(r) and the angular position 
kr in the cycle are noted. During the first half of the cycle, 
there is a shock wave on the upper surface of the airfoil 
while the flow over the lower surface remains subcritical 
throughout the entire cycle. The pressure distributions in- 
dicate that the shock position oscillates over approximately 
12% of the chord along the upper surface, requiring the 
spatial adaption procedure to track the movement and de- 


Upper surface - Calculated 

— — Lower surface - Calculated 


O Upper surface - Experiment 
□ Lower surface - Experiment 



Fig. II Comparison of instantaneous pressure distributions for the NACA 0012 airfoil pitching harmonically at Af ^ 
= 0.599, <*o = 4.86°, rv i = 2.44°, and k = 0.0814 computed using spatial adaption. 
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a(x) = -1.25°, kt = 210° a(x) = -2.4 1°, kx = 255° a(x) = -2.00°, kx = 307° a(x) = -0.54°, kx = 347° 


Fig. 12 Instantaneous meshes produced by the spatial adaption procedure for the NACA 0012 airfoil pitching 
harmonically at = 0.755, » 0 = 0.016°, 0l = 2.51°, and k = 0.0814. 



Fig. 13 Instantaneous contour lines (A p = 0.02) from the spatially adapted solution for the NACA 0012 airfoil 
pitching harmonically at M N = 0.755, o 0 = 0.016°, «, = 2.51°, and k = 0.0814. 
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velopment of the transient shock wave as it appears and 
disappears during the cycle. In general, the calculated results 
using the spatial adaption procedures compare well with the 
experimental data. 

The computational savings obtained by using the spatial 
adaption procedures can be estimated by comparing the 
total number of cells marched in time with the estimated 
number of cells marched in time if a globally fine mesh 
were used. The total number of cells marched in lime for the 
globally fine mesh is computed by multiplying the number 
of cells in the starting coarse mesh by the number of cycles, 
the number of time-steps per cycle, and four raised to the 
number of enrichment levels (since each cell is divided into 
four smaller cells). Therefore, for AGARD case 3, a total of 
276 million cells were marched compared to 6966 million 
(3628 x 3 x 10000 x 4 3 ) cells if a globally fine mesh were 
used, resulting in a computational savings of a factor of 
twenty-five. This factor does not include the computational 
overhead of the spatial adaption procedure, however, which 
was approximately 7% of the total CPU time. 

AGARD Case 5 

Results were obtained for the airfoil pitching with an 
amplitude of ■ 2.51° at A/ TO = 0.755 and o 0 = 0.016° 
(referred to as AGARD case 5). Figure 1 2 shows the instan- 
taneous adapted meshes and Fig. 13 shows the correspond- 


ing instantaneous density contour lines (A p = 0.02). The 
instantaneous meshes and density contour lines during the 
third cycle of motion were plotted at eight points in time. 
In each plot, the instantaneous pitch angle a(r) and the in- 
stantaneous angular position kr in the cycle are noted. The 
instantaneous meshes (Fig. 12) clearly indicate the enrich- 
ment in regions near the shock waves and near the stagnation 
points. They also show coarsened regions where previously 
enriched regions have relatively small flow gradients. The 
density contours during the cycle (Fig. 13) demonstrate the 
ability of the spatial adaption procedures to produce sharp 
transient shock waves. The corresponding surface pressure 
distributions during the third cycle of motion are shown in 
Fig. 14 for comparison with experimental data. 23 In each 
pressure plot the instantaneous pitch angle a{r) and the an- 
gular position Jtr in the cycle are noted. During the first part 
of the cycle there is a shock wave on the upper surface of 
the airfoil, and the flow over the lower surface is predomi- 
nately subcritical. During the latter part of the cycle the flow 
about the upper surface is subcritical while a shock forms 
along the lower surface. The pressure distributions indicate 
that the shock position oscillates over approximately 25% 
of the chord along the upper and lower surfaces, requiring 
the spatial adaption procedure to accurately track the move- 
ment and development of the transient shocks. In general 
the calculated results using the spatial adaption procedures 
compare well with the experimental data. Figure 15 shows 


— — Upper surface - Calculated 
— — Lower surface - Calculated 


O Upper surface - Experiment 
□ Lower surface - Experiment 





Fig. 14 Comparison of instantaneous pressure distributions for the NACA 0012 airfoil pitching harmonically at 
= 0.755, c*o = 0.016°, ori = 2.51°, and k = 0.0814 computed using spatial adaption. 
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Fig. 15 Variation of number of cells throughout a 
cycle of motion for the NACA 0012 airfoil 
pitching harmonically at M = 0.755, a 0 = 
0.016°, oi = 2.51°, and k = 0.0814. 

the variation of the number of cells in the mesh throughout 
the third cycle of motion. In this figure, the maximum 
number of cells is approximately 14400 and the minimum 
number of cells is near 12900. These values indicate a sig- 
nificant computational savings when spatial adaption is used 
when compared to using a globally fine mesh of 232192 
cells (comparable mesh density of 3 levels of enrichment on 
starting coarse mesh of 3628 cells is 3628 x 4 3 = 232192 
cells) for similar spatial accuracy. 

The computational savings obtained by using the spatial 
adaption procedures can again be estimated by comparing 
the total number of cells marched in time with the number of 
cells marched in time if a globally fine mesh were used. For 
AGARD case 5, a total of 413 million cells were marched 
compared to 6966 million cells for the globally fine mesh 
resulting in a computational savings of a factor of 17. This 
factor shows a decrease in savings compared to the savings 
factor obtained for AGARD case 3 due to the presence of 
a shock wave during the majority of the cycle of motion. 
Further, the computational savings factor docs not include 
the computational overhead of the spatial adaption procedure 
which was approximately 7.5% of the total CPU time. 


Concluding Remarks 

Spatial adaption procedures for the accurate and effi- 
cient solution of steady and unsteady inviscid flow problems 
were described. The adaption procedures were developed 
and implemented within a two-dimensional unstructured- 
grid upwind-type Euler code. These procedures involve 
mesh enrichment and mesh coarsening to either add points 
in high gradient regions of the flow or remove points where 
they arc not needed, respectively, to produce solutions of 
high spatial accuracy at minimal computational cost. A 


novel approach for detecting flow features based on the sub- 
stantial derivative of density was used as an enrichment in- 
dicator. This enrichment indicator worked well for detecting 
developing shock waves in unsteady flows. This is a signifi- 
cant improvement over the more commonly used enrichment 
indicators based on the instantaneous solution (such as first 
or second differences in density) that miss the initial shock 
wave formation, especially for cases where the shock waves 
periodically appear and disappear in time. 

Steady and unsteady transonic results were presented 
for the NACA 0012 airfoil to demonstrate applications of the 
spatial adaption procedures to two-dimensional problems. 
The unsteady flow results were obtained for an airfoil pitch- 
ing harmonically about the quarter chord. Both the steady 
and unsteady solutions obtained using spatial adaption were 
shown to be of high spatial accuracy, primarily in that the 
shock waves were very sharply captured. Comparing the 
cost of results obtained on a spatially adapted mesh with the 
estimated cost of results obtained on a globally fine mesh of 
comparable mesh density, a computational savings of a fac- 
tor of approximately fifty-three was achieved for the steady 
calculations using four levels of enrichment. Similar calcu- 
lations using structured grid algorithms were estimated to be 
an order of magnitude more expensive for comparable spatial 
accuracy. Comparisons of coefficients from the spatial adap- 
tion solution with those of Pulliam and Barton, 20 indicated 
agreement to within 1% of the lift and drag coefficients and 
to within 5% for the moment coefficient. Similarly, for the 
unsteady calculations, a computational savings of a factor of 
as much as twenty-five was achieved and comparisons that 
were made with experimental pressure data indicated good 
agreement. 
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